Adaptive Mesh Strategy for Efficient Use of Interface Elements in a 3D Probabilistic Explicit Cracking Model for Concrete

In this paper, the development of a 3D adaptive probabilistic explicit cracking model for concrete is reported. The contribution offered herein consists in a new adaptive mesh strategy designed to optimize the use of interface elements in probabilistic explicit cracking models. The proposed adaptive mesh procedure is markedly different from other strategies found in the literature, since it takes into account possible influences on the redistribution of stresses after cracking and can also be applied to purely deterministic cracking models. The process of obtaining the most appropriate adaptive mesh procedure involved the development and evaluation of three different adaptivity strategies. Two of these adaptivity strategies were shown to be inappropriate due to issues related to stress redistribution after cracking. The validation results demonstrate that the developed adaptive probabilistic model is capable of predicting the scale effect at a level similar to that experimentally observed, considering the tensile failure of plain concrete specimens. The results also show that different softening levels can be obtained. The proposed adaptive mesh strategy proved to be advantageous, being able to promote significant reductions in the simulation time in comparison with the classical strategy commonly used in probabilistic explicit cracking models.


Introduction
Several studies have been conducted using probabilistic numerical models that deal directly with the heterogeneity of concrete [1][2][3][4][5][6][7][8][9][10][11][12][13].With the statistical distribution of mechanical properties and possible use of fracture mechanics concepts, the failure process of structures can be analyzed in a more realistic and comprehensive way.In such models, manifestations of the scale effect, which involve the variation in global mechanical property evaluations with the change in the structure size, and softening behavior (gradual loss of load-bearing capacity) are intrinsic consequences, and the statistical treatment of global results is enabled by the Monte Carlo method.
One of the probabilistic approaches directed to the modeling of concrete behavior is the explicit cracking approach originally proposed by Rossi and Richer [1].This approach is based on the finite element method (FEM), being mainly characterized by the use of zero-thickness interface elements to which the following purposes are assigned: to receive mechanical properties distributed randomly and to represent the material discontinuities (failure planes) explicitly.The models developed from this idea have advantages in situations of crack opening assessment, which can be performed directly, in addition to favoring the implementation of strategies to deal directly with phenomena that occur during the fracture process, such as the contact between crack surfaces.
The development of 3D probabilistic explicit cracking models remains at the initial stage, with considerable limitations imposed by their high computational cost.The inclusion of interface elements in the mesh significantly increases the number of nodes and, consequently, the execution time of the analyses, especially when quadratic elements are employed.It is important to think of strategies to deal with this issue, since threedimensional analyses are indispensable in some cases, with the cracking process tending to be naturally three-dimensional.In such cases, for more realistic and comprehensive modeling, the fracture mechanisms cannot be neglected in any dimension.For safety and durability reasons, an important potential application of 3D models concerns, for example, the study of thermal cracking of massive concrete structures, such as dams and nuclear power plants.The computational cost issue can be handled by using a finite element code designed for high performance, proper data structures [14,15], and mesh adaptivity and computational parallelization strategies.
Mota et al. [11] developed a 3D probabilistic explicit cracking model capable of reproducing different levels of softening behavior and predicting the scale effect in tension with an intensity similar to that experimentally observed.This is of great importance, since the scale effect and softening behavior occur naturally for quasi-brittle materials such as concrete and strongly influence the fracture process of these materials.The mentioned model was implemented following the classical procedure of using interface elements between all the solid elements since the beginning of the analysis.Although the code was written on a high-performance FEM platform and parallelization strategies were employed, the computational cost can still be considerably reduced by adopting a mesh adaptivity strategy for the controlled use of interface elements.
As the use of interface elements is a strong time-consuming factor, an adaptivity strategy for inserting interface elements only in the regions of the mesh where they are indispensable is an important alternative to speed up the analyses.There are studies performed with adaptive mesh procedures in accordance with this idea [16][17][18][19][20][21][22][23], but these works do not encompass probabilistic cracking models and almost all of them [16][17][18][19][20][21][22] do not take into account the influence of the proposed adaptivity strategy on the redistribution of stresses after cracking.Furthermore, the adaptivity strategies proposed by Zhan and Meschke [22] and Rodrigues et al. [23] consider a different type of interface element (interface solid element with a high aspect ratio).Thus, the adaptive mesh strategies developed in all the mentioned studies [16][17][18][19][20][21][22][23] cannot be applied to the probabilistic explicit cracking approach, so a new strategy is necessary.An adaptivity procedure appropriate for probabilistic explicit cracking models can prove to be very efficient, for example, in simulations in which the structural failure occurs due to concentrated cracking in delimited regions of the mesh.
In view of the foregoing considerations, this work reports the development of a 3D adaptive probabilistic cracking model based on the explicit cracking approach proposed by Rossi and Richer [1].The model in question represents an improvement of the one implemented by Mota et al. [11].The contribution offered herein consists in an adaptive mesh strategy designed to optimize the use of interface elements in probabilistic explicit cracking models, taking into account possible influences on the redistribution of stresses after cracking.This strategy is markedly different from the strategies found in the literature [16][17][18][19][20][21][22][23] and can also be applied to purely deterministic cracking models that employ interface elements.

Empirical Relations for Scale Effect
Rossi et al. [24], by means of an extensive experimental investigation, obtained analytical expressions that define correlations between concrete heterogeneity and scale effects, contributing to the validation of probabilistic models.These expressions depend on the V T /V A ratio, which is a parameter of heterogeneity, and on the compressive strength ( f c ).The parameter V T is the total volume of the specimen, and V A is the volume of the coarsest aggregate, modeled by a sphere.The expressions are applicable to any concrete, except fiber-reinforced and lightweight concretes, with an f c ranging from 35 MPa to 130 MPa and with the V T /V A ratio varying from 10 to 10,000 in order of magnitude.For the mean value of tensile strength or mean tensile strength ( f mean t ) and the coefficient of variation (CV), which is a measure of dispersion defined as the ratio between the standard deviation (SD) and the mean tensile strength, Rossi et al. [24] obtained the following relations: (1) in which a = 6.5 MPa, A = 0.35, and the values of γ and β are defined by γ = 0.25 − 3.6 × 10 −3 f c c 1 with c 1 = 1 MPa. Figure 1 portrays the diagrams of γ and β versus f c .With the aid of these diagrams, it is possible to note that, with γ and β being always positive, the higher the V T /V A ratio in Equations ( 1) and (2), the lower the mean tensile strength and the coefficient of variation.Figure 1 also shows that, with the compressive strength being a quality indicator of the cement paste, when f c increases, γ decreases, leading to higher values of mean tensile strength.On the other hand, higher values of f c result in higher values for β and, consequently, in lower values for CV.strength ( ).The parameter  is the total volume of the specimen, and  is the volume of the coarsest aggregate, modeled by a sphere.The expressions are applicable to any concrete, except fiber-reinforced and lightweight concretes, with an  ranging from 35 MPa to 130 MPa and with the   ⁄ ratio varying from 10 to 10,000 in order of magnitude.For the mean value of tensile strength or mean tensile strength ( ) and the coefficient of variation ( ), which is a measure of dispersion defined as the ratio between the standard deviation ( ) and the mean tensile strength, Rossi et al. [24] obtained the following relations: with  = 1 MPa. Figure 1 portrays the diagrams of  and  versus  .With the aid of these diagrams, it is possible to note that, with  and  being always positive, the higher the   ⁄ ratio in Equations ( 1) and (2), the lower the mean tensile strength and the coefficient of variation.Figure 1 also shows that, with the compressive strength being a quality indicator of the cement paste, when  increases,  decreases, leading to higher values of mean tensile strength.On the other hand, higher values of  result in higher values for  and, consequently, in lower values for .1)-( 4)).
According to Rossi et al. [24], the tensile strength dispersion is occasioned by the material heterogeneity, being directly associated with the scale effect.Small values of , obtained either by increasing the   ⁄ ratio or by increasing the value of  , represent a low influence of the heterogeneity on the results and, consequently, a low-intensity scale effect.

3D Adaptive Probabilistic Explicit Cracking Model
Figure 2 outlines the 3D probabilistic explicit cracking model without mesh adaptivity developed by Mota et al. [11], which is the basis of the present work, including for comparison purposes.The main characteristics of this model are described as follows: According to Rossi et al. [24], the tensile strength dispersion is occasioned by the material heterogeneity, being directly associated with the scale effect.Small values of CV, obtained either by increasing the V T /V A ratio or by increasing the value of f c , represent a low influence of the heterogeneity on the results and, consequently, a low-intensity scale effect.

3D Adaptive Probabilistic Explicit Cracking Model
Figure 2 outlines the 3D probabilistic explicit cracking model without mesh adaptivity developed by Mota et al. [11], which is the basis of the present work, including for comparison purposes.The main characteristics of this model are described as follows: • The three-dimensional treatment is basically given by the use of tetrahedral solid elements of 10 nodes; • Zero-thickness interface elements of 6 + 6 nodes are previously inserted between all the solid elements in accordance with the classical method of mesh formation; • Quadratic interpolation functions are used in the elements to obtain results with a higher precision; • The heterogeneous and probabilistic character is established by the random assignment of tensile strength values to the interface elements according to the Weibull distribution.Thus, concrete is being defined as a heterogeneous material, since there is no tensile strength variation for a homogenous material.Weibull's statistical law was chosen for two reasons: (1) it was developed for brittle materials, which can be associated with the fact that a brittle law is considered for the material locally; (2) negative values of random properties are not possible; • Solid elements have a linear elastic behavior;

•
Interface elements follow a sufficiently rigid and brittle behavior under tension, which is naturally associated with the Rankine failure criterion, resulting in local geometric and material discontinuity when the tensile strength of a given interface element is reached, causing the occurrence of a cracked interface element.Thus, during the computational analysis, there are local nonlinearities that generate nonlinearity in the global response; • Structural softening is naturally obtained according to the distribution of tensile strength among the interface elements; • When there is contact between the faces of cracked interface elements (contact problem locally), there is a reactivation of stiffness in the interface elements to take into account compression and friction effects; • The acquisition of structural responses occurs by means of the Monte Carlo method.

•
The three-dimensional treatment is basically given by the use of tetrahedral solid elements of 10 nodes; • Zero-thickness interface elements of 6 + 6 nodes are previously inserted between all the solid elements in accordance with the classical method of mesh formation; • Quadratic interpolation functions are used in the elements to obtain results with a higher precision;

•
The heterogeneous and probabilistic character is established by the random assignment of tensile strength values to the interface elements according to the Weibull distribution.Thus, concrete is being defined as a heterogeneous material, since there is no tensile strength variation for a homogenous material.Weibull's statistical law was chosen for two reasons: (1) it was developed for brittle materials, which can be associated with the fact that a brittle law is considered for the material locally; (2) negative values of random properties are not possible; • Solid elements have a linear elastic behavior; • Interface elements follow a sufficiently rigid and brittle behavior under tension, which is naturally associated with the Rankine failure criterion, resulting in local geometric and material discontinuity when the tensile strength of a given interface element is reached, causing the occurrence of a cracked interface element.Thus, during the computational analysis, there are local nonlinearities that generate nonlinearity in the global response; • Structural softening is naturally obtained according to the distribution of tensile strength among the interface elements; When there is contact between the faces of cracked interface elements (contact problem locally), there is a reactivation of stiffness in the interface elements to take into account compression and friction effects; The acquisition of structural responses occurs by means of the Monte Carlo method.The constitutive relation for the interface element can be basically written in the form of the following damage-type model [11]: in which δ is the damage parameter, k 11 and k 22 are penalty parameters associated with the tangential stiffness, k 33 is related to the normal stiffness, w t1 and w t2 are tangential relative displacements, and w n is the normal relative displacement.In order to obtain a brittle behavior, it is assumed that δ = 0 for the intact interface element, and δ = 1 when the tensile stress reaches the tensile strength (occurrence of a cracked interface element).
The friction between crack surfaces, which are represented by faces of cracked interface elements, is taken into account by means of the Mohr-Coulomb criterion [25] without cohesion: where τ is the shear stress or, more specifically, the acting friction stress in the cracked interface element; σ c is the absolute value of the normal compressive stress; and ϕ is the angle of friction.A more detailed description of the friction model is provided by Mota et al. [12].The next subsections describe the implementation of mesh adaptivity in the model portrayed in Figure 2.For convenience, the regions where interface elements can be inserted are referred to as interfacial regions of solid elements or simply interfacial regions.

Generation and Distribution of Random Tensile Strength
For the Weibull distribution [26], the cumulative distribution function and the inverse cumulative distribution function related to a random variable x can be written respectively as G(y, b, c) = c(−ln(1 − y)) in which x ≥ 0, with this variable corresponding to the tensile strength in the present paper; b is the shape parameter, associated with the coefficient of variation; and c is the scale parameter, directly related to the mean value of x.
The generation of random values of tensile strength is based on the inverse transformation method or uniform probability transformation technique [27].These values are generated and distributed to interfacial regions of solid elements according to Algorithm A1 (Appendix A).If the values of b and c are fixed, a different Monte Carlo sample for the same material will be obtained each time Algorithm A1 is used.

Adaptive Mesh Strategies for Efficient Use of Interface Elements
Regarding finite element meshes for mechanical problems, a relevant challenge is how to represent the crack explicitly.For a mesh composed only of solid elements, a crack at an interfacial region of the elements, as indicated in Figure 3a, can be detected relatively simply through the calculation of nodal stresses in the region.However, the incorporation of the geometric discontinuity generated by the crack requires strategies of greater complexity.
Figure 3b illustrates the classical strategy adopted in probabilistic explicit cracking models to deal with the appearance of cracks, while Figure 3c-e present the three adaptive mesh strategies developed and evaluated in the context of this work.Figure 4 provides more details about the adaptivity options, with the solid elements being displayed in a reduced size in the crack zone so that the interface elements are highlighted.Two dimensional representations are used for simplification purposes.
the mesh is initially prepared with interface elements inserted between all the solid elements.Thus, the crack observed in Figure 3a is detected directly in the interface element, which produces a geometric discontinuity when its rigidity is nullified, instantly generating stress redistribution around the crack.However, this strategy causes a significant increase in the number of nodes from the beginning of the analysis, with interface elements being used even in regions where these elements may be unnecessary.In the here-called Classical Strategy of Mesh Formation (CSMF), shown in Figure 3b, the mesh is initially prepared with interface elements inserted between all the solid elements.Thus, the crack observed in Figure 3a is detected directly in the interface element, which produces a geometric discontinuity when its rigidity is nullified, instantly generating stress redistribution around the crack.However, this strategy causes a significant increase in the number of nodes from the beginning of the analysis, with interface elements being used even in regions where these elements may be unnecessary.
The three previously mentioned mesh adaptivity strategies were developed so that interface elements could be employed only in regions where they are indispensable.In these strategies, the mesh is initially composed only of solid elements, with interface elements being used during the analysis as cracking occurs.Consequently, after the beginning of the analysis, there may or may not be interface elements in the interfacial regions of the solid elements.When a crack is detected at an interfacial region with no interface element, a cracked interface element is inserted in this region.However, for mesh conformity, it is also necessary to place intact interface elements in the vicinity of the cracked interface element.In the first mesh adaptivity strategy, referred to as Adaptive Mesh Strategy 1 (AMS1) and portrayed in Figures 3c and 4a, intact interface elements are inserted around the solid element that provides the top face of the cracked interface element.As the choice for the solid element related to the top face was purely arbitrary, one could opt to use intact interface elements around the solid element associated with the base face of the cracked interface element instead.In a 2D mesh composed of triangular elements, this strategy would entail the inclusion of at most two intact interface elements for each detected crack.In a 3D configuration with tetrahedral solid elements, at most three intact interface elements are included for each new crack.
In Adaptive Mesh Strategy 2 (AMS2), depicted in Figures 3d and 4b, intact interface elements are placed around the solid elements that provide the base and top faces of the cracked interface element.In this way, for each new crack, a maximum of four intact interface elements would be inserted in a 2D case with triangular elements; on the other hand, in 3D meshes of tetrahedrons, a maximum of six intact interface elements can be added.
In Adaptive Mesh Strategy 3 (AMS3), which is represented in Figures 3e and 4c, intact interface elements are placed around all the solid elements in the vicinity of the cracked interface element.In this case, for each new crack, the maximum number of intact interface elements that can be inserted depends on the level of mesh discretization in the crack zone.
Details of the computational implementation of AMS1, AMS2 and AMS3, such as utilized arrays, explained synthesized algorithms, and further information about the insertion of cracked and intact interface elements in the mesh, are described in Appendix B.

Computational Program for Finite Element Analyses
The finite element code, written in Fortran language, was implemented in a highperformance computational platform developed by Ribeiro and Ferreira [15].This platform allows for the use of algorithms and data structures designed for high performance [14,15], having already been used as a basis for several works [28][29][30][31].Details related to the developed computational program are described in Appendix C, with an overview of the implemented code being presented in Algorithm A5.This appendix includes the description of a method adopted to favor proper stress redistributions during the cracking process and to avoid any dependence of results on the load increment size, called the most dangerous element method.
The application of the Monte Carlo method requires Algorithm A5 to be executed a number of times, corresponding to a number of Monte Carlo samples whose loaddisplacement responses may naturally vary considerably from one sample to another.In order to deal automatically with this variation, two simulation stopping criteria were implemented, taking into account some load-displacement diagram information: maximum load (F max ), ultimate load (F u ), and ultimate displacement (δ tot ), as indicated in Figure 5.The simulation stopping criteria are named as follows: SSC1, whose evaluation parameter is the ratio between  and  , with  being necessarily less than  ; and SSC2, whose evaluation parameter is simply  .It is also possible to opt for both criteria simultaneously.More specifically, these criteria are defined as follows: The simulation stopping criteria are named as follows: SSC1, whose evaluation parameter is the ratio between F u and F max , with F u being necessarily less than F max ; and SSC2, whose evaluation parameter is simply δ tot .It is also possible to opt for both criteria simultaneously.More specifically, these criteria are defined as follows: • SSC1 [ψ a ] → the criterion is met, and the simulation is terminated, when → the criterion is met, and the simulation is terminated, when δ tot ≥ ψ b .with the limit parameters ψ a and ψ b being set according to the particularities of the considered analysis, so that 0 ≤ ψ a < 1 and ψ b > 0.
It is important to point out that values of ψ a close to 1 should be avoided, since such values do not guarantee the effective knowledge of the peak load.As load fluctuations can naturally occur in the load-displacement diagram, inappropriate values of ψ a can easily cause false peak-load results.For each type of mechanical problem, it is necessary to carry out previous analyses to determine suitable values of ψ a .

Analyzed Cases
The numerical simulations were divided into four stages: (1) mechanical evaluation of the adaptive mesh strategies; (2) inverse analysis; (3) validation; (4) assessment of time savings.Figures 6-8 show the geometries and meshes of the analyzed cases, with the details presented in Table 1.In this table, the columns related to interface elements and nodes inform the maximum possible number of these entities in the mesh, remembering that, in case of adaptive mesh use, the informed numbers are not necessarily reached at the end of the analysis; the parameter 2V e denotes the average value of the sum of the volumes of the two solid elements that provide the base and the top faces of the interface elements.In stage 1, direct tensile test simulations were performed by using the square-based prisms (P1, P2, and P3) shown in Figure 6, with P3 corresponding to a cube.Boundary and loading conditions in agreement with the simulated test have been considered, with displacement restrictions being imposed on the base nodes, and displacement increments being applied on the top nodes in the longitudinal direction.Mechanical responses obtained by using the adaptive mesh strategies were compared with the responses provided by the classical strategy of mesh formation.For P1 and P2, the comparisons refer to values of the average, standard deviation, and coefficient of variation of tensile strength,    In stage 1, direct tensile test simulations were performed by using the square-based prisms (P1, P2, and P3) shown in Figure 6, with P3 corresponding to a cube.Boundary and loading conditions in agreement with the simulated test have been considered, with displacement restrictions being imposed on the base nodes, and displacement increments being applied on the top nodes in the longitudinal direction.Mechanical responses obtained by using the adaptive mesh strategies were compared with the responses provided by the classical strategy of mesh formation.For P1 and P2, the comparisons refer to values of the average, standard deviation, and coefficient of variation of tensile strength, In stage 1, direct tensile test simulations were performed by using the square-based prisms (P1, P2, and P3) shown in Figure 6, with P3 corresponding to a cube.Boundary and loading conditions in agreement with the simulated test have been considered, with displacement restrictions being imposed on the base nodes, and displacement increments being applied on the top nodes in the longitudinal direction.Mechanical responses obtained by using the adaptive mesh strategies were compared with the responses provided by the classical strategy of mesh formation.For P1 and P2, the comparisons refer to values of the average, standard deviation, and coefficient of variation of tensile strength, and to the global load-displacement responses.For P3, the comparisons concern the local results of stress redistribution after crack appearance.The purpose of these comparisons was to identify the most suitable adaptive mesh strategy to be adopted in stages 2, 3, and 4, using the classical strategy of mesh formation (CSMF) to provide reference results.The mean tensile strength of each set of Monte Carlo samples ( f mean t ) and the corresponding standard deviation (SD) are calculated by the following expressions [27]: (10) in which M is the number of Monte Carlo samples, and f max is the maximum load of the i-th Monte Carlo sample, and A is the cross-sectional area.
In stages 2 and 3, the simulations with P1 and P2 were based on the experimental investigation conducted by Rossi et al. [24], which was chosen for the following reasons: • The study refers to concrete behavior in a direct tensile test, in which the transition from diffuse to localized cracking is the most important phenomenon and the most difficult to capture; • Concretes with different mechanical properties were produced, enabling the evaluation of different scenarios; • A large number of experimental samples were tested, implying statistical relevance; • Experimental results are expressed using Equations ( 1) and (2) with great accuracy, in such a way that it is completely possible to use these equations to obtain expected values according to the scale effect.
In stage 2, inverse analysis procedures described by Mota et al. [11] were carried out with the purpose of determining the parameters b and c (statistical parameters related to Equations ( 7) and ( 8)) for P1 and B1.It is important to highlight that these parameters naturally depend on the volume of the mesh elements, in a manner similar to that in which the values of f mean t and CV given by Equations ( 1) and (2) depend on the volume of the considered specimens.The parameters b and c are scale effect parameters at the local level, while the values of f mean t and CV provided by Equations ( 1) and ( 2) are global-level scale effect results (specimen responses).
Basically, each inverse analysis procedure corresponds to analyzing the same case with different pairs of b and c to define the most appropriate pair of such parameters, taking into account empirical results.A considerable number of Monte Carlo processes is necessary.One Monte Carlo process provides, for example, the mean and the standard deviation of a certain global mechanical response for a given pair of b and c.Repeating this procedure for different pairs of b and c, two surfaces are obtained: one for the mean and another for the standard deviation.These surfaces are properly fitted using polynomial expressions.The expected results for the mean and standard deviation, which can be provided using Equations ( 1) and (2) for simulations with the prism P1, can be considered as horizontal planes on the b-c domain.In the next step, the intersection curves between these planes and the fitted surfaces must be determined.Figure 9 shows these curves, considering that the tensile strength of the samples, which is given by the maximum load divided by the cross-sectional area, is the mechanical response to be observed.In this figure, the intersection point of the curves represents the most suitable pair of b and c.More details related to the inverse analysis procedures are presented by Mota et al. [11].For beam B1, detailed in Figure 7 and experimentally tested by Hordijk [32], the boundary and loading conditions were consistent with the simulated four-point bending test, with points A and B indicating the locations of displacement restrictions and points C and D signaling the places of application of displacement increments.Point E marks the deflection measurement location.The acting load was measured as the sum of the forces applied to the supports that contain points C and D. The possibility of using interface elements was restricted to the region highlighted in red in Figure 8.In the simulated test, due to its particularities, the failure is due to tension.
In stage 3, the ability to predict the scale effect was investigated.The values of  and  obtained from an inverse analysis with P1 were used in simulations with P2, which has larger dimensions, for the same concretes.In addition to depending on the material, the parameters  and , as previously mentioned, depend on the volume of the mesh solid elements.Therefore, the meshes for P1 and P2 were prepared keeping the same average size of elements.In this way, it was possible to evaluate only the effect of increasing the size of the simulated specimen.
Also due to the dependency between the parameters  and  and the size of the elements, the mesh for B1 was generated with the elements of the central region having a very low level of volumetric variation.In relation to P3, there was no need to guarantee For beam B1, detailed in Figure 7 and experimentally tested by Hordijk [32], the boundary and loading conditions were consistent with the simulated four-point bending test, with points A and B indicating the locations of displacement restrictions and points C and D signaling the places of application of displacement increments.Point E marks the deflection measurement location.The acting load was measured as the sum of the forces applied to the supports that contain points C and D. The possibility of using interface elements was restricted to the region highlighted in red in Figure 8.In the simulated test, due to its particularities, the failure is due to tension.
In stage 3, the ability to predict the scale effect was investigated.The values of b and c obtained from an inverse analysis with P1 were used in simulations with P2, which has larger dimensions, for the same concretes.In addition to depending on the material, the parameters b and c, as previously mentioned, depend on the volume of the mesh solid elements.Therefore, the meshes for P1 and P2 were prepared keeping the same average size of elements.In this way, it was possible to evaluate only the effect of increasing the size of the simulated specimen.
Also due to the dependency between the parameters b and c and the size of the elements, the mesh for B1 was generated with the elements of the central region having a very low level of volumetric variation.In relation to P3, there was no need to guarantee low volumetric variation among the elements, since there was no inverse analysis for this case.For this reason, the value of 2V e is not important information for P3, as indicated in Table 1.
In stage 4, the efficiency of the adaptive mesh strategy selected in stage 1 was investigated.In order to do this, the simulation times for the analyses performed with the adaptive mesh strategy and the CSMF were compared.As shown in Table 1, the results for P1, P2, and B1 were taken into account in this stage.
Table 2 displays the compressive strength, the modulus of elasticity (E), and the volume of the coarsest aggregate of the considered concretes, as well as the heterogeneity parameters V T /V A and 2V e /V A .Concretes C1, C2, and C3, experimentally studied by Rossi et al. [24], were applied to P1 and P2.Concrete C4, related to the experimental investigation conducted by Hordijk [32], was applied to B1, for which the value of 2V e /V A corresponds to the region in red in Figure 8. Table 3 contains the values of complementary modeling parameters: the angle of friction, which corresponds to an average estimate based on studies found in the literature [33,34]; the Poisson's ratio of concretes C1, C2, C3, and C4; and the Poisson's ratios and moduli of elasticity for P3 and for the steel support of B1.The adopted value for the penalty parameters k 11 , k 22 , and k 33 related to the constitutive relation of intact interface elements was 1 × 10 10 MN/dm 3 .This value was chosen from a calibration for which the results of a numerical linear analysis were similar to the analytical ones, without causing problems in the solution.The nomenclature applied to the Monte Carlo samples involves the geometry, mesh, and concrete in question, according to the terminology used in Tables 1 and 2. Altogether, there are seven classes of Monte Carlo samples: P1-C1, P1-C2, P1-C3, P2-C1, P2-C2, P2-C3, and B1-C4.
It is worth mentioning that, for the considered probabilistic cracking model, there are some restrictions for performing a mesh sensitivity analysis in which a convergence of results is evaluated as the degree of mesh refinement increases.The variation in the degree of mesh refinement naturally implies a change in the volume of the solid elements.However, the values of b and c depend on the volume of the solid elements (as mentioned before), which means that different values of b and c need to be used when the volume of the solid elements changes, modifying the characteristics of the problem and compromising the mesh sensitivity analysis.A mesh sensitivity analysis can be performed considering the elastic response (in order to ensure the correct elastic response, which is guaranteed in the cases analyzed in this paper) or considering the problem without tensile strength variation (with the same value of tensile strength for all interface elements in the model).In the present paper, the degree of mesh refinement is not a problem, since the adopted meshes provide a level of cracking information that enables the investigation described herein to be carried out.If results with a very high level of cracking information were needed, the meshes could be more refined (which would naturally generate a higher computational cost).But this is not the case in this work, whose purpose is to present a new adaptive mesh strategy and its viability, which does not require results with a very high level of cracking information.

Results
The results of the analyses conducted in the stages described above are exposed and discussed in the next subsections.

Mechanical Evaluation of the Adaptive Mesh Strategies
For the evaluation of the mesh adaptivity strategies, a total of 80 Monte Carlo samples were simulated for each adaptive mesh strategy as well as for the CSMF.Three pairs of values of b and c that represent different material characteristics were chosen based on previous analyses carried out only to obtain coherent initial values of these parameters: b = 1.2 and c = 12.0 MPa; b = 3.0 and c = 8.0 MPa; and b = 20.0 and c = 7.5 MPa.The following parameters were calculated: mean tensile strength, obtained via Equation ( 9); standard deviation, determined using Equation (10); coefficient of variation (CV), given by the ratio between SD and f mean t ; and the discrepancy between f mean t associated with the adaptive mesh strategy in question and f mean t obtained by using the CSMF (∆ 1 ).The results for P1 and P2 are presented in Tables 4 and 5, respectively.Tables 4 and 5 show that, for the adaptive mesh strategies, the value of f mean t decreases as the number of interface elements inserted for each new crack increases.In other words, the highest and the lowest value of f mean t are generated with AMS1 and AMS3, respectively.This was observed for the six classes of samples considered in the tables in question.In relation to the comparison with the results of f mean t obtained by using the CSMF, for samples related to P1, Table 4 reveals that the discrepancy ∆ 1 was between 16.4% and 57.5% for AMS1, between 14.5% and 34.3% for AMS2, and between −0.2% and 1.1% for AMS3.For samples related to P2, Table 5 shows that ∆ 1 ranged between 21.7% and 58.3% for AMS1, between 19.9% and 37.0% for AMS2, and between −0.1 and 2.2% for AMS3.
The results contained in Tables 4 and 5 indicate that only AMS3 is capable of generating values of f mean t similar to those produced with the CSMF.Considering the six classes of simulated samples, the maximum absolute value of ∆ 1 for AMS3, equal to 2.2%, is significantly lower than the value of ∆ 1 referring to AMS1 and AMS2.Furthermore, the values of SD and CV for AMS3 also have a pronounced similarity with those associated with the CSMF.It is also important to highlight the occurrence of the scale effect on the values of f mean t , SD, and CV for all the considered strategies.In regard to f mean t , it can be seen in Tables 4 and 5 that, for AMS3 for example, the value of this parameter changed from 2.95 MPa in P1-C1 to 2.85 MPa in P2-C1, from 3.77 MPa in P1-C2 to 3.57 MPa in P2-C2, and from 6.13 MPa in P1-C3 to 5.87 MPa in P2-C3.However, scale effect occurrences will be rigorously evaluated in Section 5.3 with values of b and c obtained from inverse analyses, which are necessary for an effective comparison with empirical results.
In order to promote a more specific assessment of the mechanical responses obtained with the adaptive mesh strategies, Figure 10 displays load-displacement diagrams for a certain Monte Carlo sample of classes P2-C1, P2-C2, and P2-C3.It should be mentioned that, for each class, the curves refer to exactly the same sample; that is, the random values of tensile strength at the interfacial regions of solid elements were kept the same for all the strategies.
Figure 10 shows that AMS3 tends to produce load-displacement responses very similar to those generated by the CSMF.As seen in Figure 10b, occasional differences can occur mainly in the post-peak branch.The reason for this is that the stress calculation referring to the interfacial regions of solid elements does not happen in the exact same way for the two strategies (see the process of obtaining the array Savg in Appendix B), which may cause different crack propagation paths depending on the random distribution of tensile strength values in the mesh.However, comparing the mechanical responses of other samples, in general, there were no accentuated differences in the post-peak branch, indicating that the CSMF and AMS3 indeed tend to generate similar fracture processes.
Regarding AMS1 and AMS2, Figure 10 demonstrates that both strategies produce loaddisplacement diagrams markedly different from those generated by the CSMF.Although there is a certain degree of similarity in the shape of the curves, the cracking localization process starts at a higher load level when AMS1 and AMS2 are employed.Consequently, the maximum loads are naturally also higher.This divergent behavior occurs because AMS1 and AMS2 do not generate an amount of intact interface elements in regions close to cracks capable of enabling a stress redistribution process similar to that observed with the CSMF, as discussed next.

Stress Redistribution Verification
In order to clarify the divergences observed in the mechanical responses obtained with AMS1 and AMS2 in relation to the results provided by the CSMF and AMS3, simulations with P3 were performed taking into account the occurrence of a previously directed crack.The redistribution of stresses in the vicinity of the crack was investigated for each strategy for interface element usage.Figure 11  Tables 4 and 5 that, for AMS3 for example, the value of this parameter changed from 2.95 MPa in P1-C1 to 2.85 MPa in P2-C1, from 3.77 MPa in P1-C2 to 3.57 MPa in P2-C2, and from 6.13 MPa in P1-C3 to 5.87 MPa in P2-C3.However, scale effect occurrences will be rigorously evaluated in Section 5.3 with values of  and  obtained from inverse analyses, which are necessary for an effective comparison with empirical results.
In order to promote a more specific assessment of the mechanical responses obtained with the adaptive mesh strategies, Figure 10 displays load-displacement diagrams for a certain Monte Carlo sample of classes P2-C1, P2-C2, and P2-C3.It should be mentioned that, for each class, the curves refer to exactly the same sample; that is, the random values of tensile strength at the interfacial regions of solid elements were kept the same for all the strategies.Figure 10 shows that AMS3 tends to produce load-displacement responses very similar to those generated by the CSMF.As seen in Figure 10b, occasional differences can occur mainly in the post-peak branch.The reason for this is that the stress calculation referring to the interfacial regions of solid elements does not happen in the exact same way for the two strategies (see the process of obtaining the array Savg in Appendix B), which may cause different crack propagation paths depending on the random distribution of tensile strength values in the mesh.However, comparing the mechanical responses of other samples, in general, there were no accentuated differences in the post-peak branch, indicating that the CSMF and AMS3 indeed tend to generate similar fracture processes.Regarding AMS1 and AMS2, Figure 10 demonstrates that both strategies produce load-displacement diagrams markedly different from those generated by the CSMF.Although there is a certain degree of similarity in the shape of the curves, the cracking localization process starts at a higher load level when AMS1 and AMS2 are employed.Consequently, the maximum loads are naturally also higher.This divergent behavior occurs because AMS1 and AMS2 do not generate an amount of intact interface elements in regions close to cracks capable of enabling a stress redistribution process similar to that observed with the CSMF, as discussed next.

Stress Redistribution Verification
In order to clarify the divergences observed in the mechanical responses obtained with AMS1 and AMS2 in relation to the results provided by the CSMF and AMS3, simulations with P3 were performed taking into account the occurrence of a previously directed crack.The redistribution of stresses in the vicinity of the crack was investigated for each strategy for interface element usage.Figure 11   After imposing a displacement of 6 × 10 −5 dm on the top face of the cube, a crack was induced in region 565.The normal stress in regions 951 and 1318 was evaluated before and right after the appearance of the crack.The results for each strategy of interface element usage are revealed in Table 6.After imposing a displacement of 6 × 10 −5 dm on the top face of the cube, a crack was induced in region 565.The normal stress in regions 951 and 1318 was evaluated before and right after the appearance of the crack.The results for each strategy of interface element usage are revealed in Table 6.It can be seen in Table 6 that the CSMF and AMS3 caused the same level of stress variation in regions 951 (from 1.2 MPa to 1.6 MPa) and 1318 (from 0.8 MPa to 0.5 MPa).This resemblance in the stress redistribution process explains why these strategies tend to produce similar mechanical results, as observed in Tables 4 and 5 and in Figure 10.
With regard to the analyses carried out with AMS1 and AMS2, Table 6 shows that region 951, without an interface element, did not undergo stress variation, with the stress before and right after the appearance of the crack remaining at a level of 1.2 MPa.On the other hand, in region 1318, which has an interface element, a considerable stress variation was noticed.
In region 951, for AMS1 and AMS2, there was no stress variation, because the change in the strain level of the solid elements in the vicinity of the crack was extremely low.It was observed that the existence of an interface element at an interfacial region close to the crack is necessary so that the stress redistribution is captured with greater sensitivity in this region.Therefore, in the analyses with the CSMF and AMS3, noticeable stress variations, captured by interface elements, were observed in the two examined regions.

Decision on Adaptive Mesh Procedure
It is assumed that the CSMF promotes adequate stress redistribution, since interface elements are previously inserted in all the interfacial regions of solid elements in this strategy, which makes the interface elements capable of capturing post-crack stress redistribution in any part of the mesh.For this reason, an appropriate and coherent adaptive mesh strategy for inserting interface elements efficiently during the analysis must provide results similar to those obtained with the CSMF.According to the evaluation of the adaptive mesh strategies, AMS1 and AMS2 cannot be used to replace the CSMF due to the occurrence of considerable divergences in the stress redistribution processes, resulting in important differences in the mechanical responses.Nevertheless, these results are of great relevance, since they can prevent these strategies from being used by other researchers in future investigations.On the other hand, the results generated with AMS3 show a high degree of equivalence with those obtained through the CSMF, which is justified by the similarity in the redistribution of stresses.Thus, considering the need to submit the mesh adaptivity to other assessments, AMS3 was chosen as the most suitable adaptive mesh strategy to continue the investigation.
AMS3 is markedly different from other adaptive mesh strategies involving interface elements found in the literature [16][17][18][19][20][21].Such a difference is mainly due to the fact that these works do not take into account the influence of the adaptive mesh strategy on the redistribution of stresses after cracking.In these works, there is no comparison between the results generated with the proposed adaptive mesh strategy and the results obtained with interface elements in all the interfacial regions of solid elements from the beginning of the analysis.
It is important to highlight that, as the probabilistic nature of the cracking model used here did not affect the development of the adaptivity strategies, AMS3 can also be applied to purely deterministic cracking models that employ interface elements.

Inverse Analysis
The values of b and c obtained through inverse analysis procedures for samples related to P1 and B1 are shown in Table 7.For the samples related to P1, values of f mean t and CV given by Equations ( 1) and (2) were taken into account.For B1-C4, the inverse analysis was based on values of the maximum load and work of fracture extracted from an experimental curve.It can be noticed that the values of b and c determined with AMS3 are very similar to those obtained using the CSMF, with an absolute discrepancy of less than 2%.Naturally, the statistical parameters are not exactly the same for the two strategies due to the possibility of small differences in the mechanical responses, as previously discussed.In order to verify the accuracy of the performed inverse analyses, the values of b and c presented in Table 7 were used to simulate 80 Monte Carlo samples for each of classes P1-C1, P1-C2, and P1-C3 and 20 Monte Carlo samples of class B1-C4.The results are displayed in Tables 8 and 9.The variable ∆ 2 in Table 8 denotes the discrepancy between the numerical modeling result referring to f mean t , obtained through Equation ( 9), and the empirical value of f mean t , determined using Equation (1).In Table 9, F exp max and W exp f represent, respectively, the experimental values of the maximum load and work of fracture until the deflection of 0.3 mm, F mean max and W mean f denote the numerical modeling results for the mean maximum load and mean work of fracture until the mentioned deflection respectively, and ∆ 3 refers to the discrepancy between F mean max and F exp max .It can be seen in both tables that all the numerical modeling results closely resemble the empirical or experimental data.The absolute values of ∆ 2 and ∆ 3 are less than 2%, associating the inverse analyses with a high degree of accuracy.Figure 12 shows the load-deflection diagrams of the 20 Monte Carlo samples (MCS) of class B1-C4, considered in Table 9, for simulations with the CSMF and AMS3.The experimental response is also displayed.As expected, although the values of b and c for B1-C4 were obtained with AMS3, these values can be normally used in simulations with the CSMF, with the mechanical responses for both strategies being quite similar, mainly in terms of the average results, as indicated in Table 9.

Validation
The validation analyses consisted in simulating 80 Monte Carlo samples for each of the classes P2-C1, P2-C2, and P2-C3, with the values of  and  presented in Table 7.The results are revealed in Table 10, with ∆ being the discrepancy between numerical modeling and empirical values for  , as previously defined.Again, a strong similarity is observed between the results provided by the CSMF and AMS3.By comparing the numerical modeling results shown in Table 10 with those displayed in Table 8, it is possible to notice the scale effect on the values of  , , and  for the three analyzed classes, which was already expected, taking into account the results discussed in Section 5.1.In Table 10, the values of ∆ range from −4.3% to 8.6%, which can be judged as a good result, considering the high complexity of the probabilistic scale effect linked to the value of  .The numerical modeling results for  and  also show good similarity with the empirical ones.Thus, it can be seen that the adaptive probabilistic model, as well as the classical alternative, is capable of predicting the scale effect at a level

Validation
The validation analyses consisted in simulating 80 Monte Carlo samples for each of the classes P2-C1, P2-C2, and P2-C3, with the values of b and c presented in Table 7.The results are revealed in Table 10, with ∆ 2 being the discrepancy between numerical modeling and empirical values for f mean t , as previously defined.Again, a strong similarity is observed between the results provided by the CSMF and AMS3.By comparing the numerical modeling results shown in Table 10 with those displayed in Table 8, it is possible to notice the scale effect on the values of f mean t , SD, and CV for the three analyzed classes, which was already expected, taking into account the results discussed in Section 5.1.In Table 10, the values of ∆ 2 range from −4.3% to 8.6%, which can be judged as a good result, considering the high complexity of the probabilistic scale effect linked to the value of f mean t .The numerical modeling results for SD and CV also show good similarity with the empirical ones.Thus, it can be seen that the adaptive probabilistic model, as well as the classical alternative, is capable of predicting the scale effect at a level similar to that experimentally observed.
For a general assessment of the mechanical responses produced with the CSMF and AMS3, Figure 13 presents the 80 load-displacement diagrams obtained for each class of Monte Carlo samples indicated in Table 10 for both strategies.With the purpose of enabling detailed examinations, Figure 14 displays the typical curves for each class until the displacement of 1.5 × 10 −3 dm.In this paper, a typical curve is one of the curves of the set of Monte Carlo sample curves that has characteristics that makes it close to an average curve.In the mentioned figures, the resemblance between the responses produced with the CSMF and AMS3 is clear, and occasional divergences between individual behaviors are easily compensated by the application of the Monte Carlo method.It is important to observe that parameter b exerts a strong influence on the results.Low values of this parameter are associated with high variability in the obtained diagrams.Moreover, the higher the value of b, the lower the level of nonlinearity in the pre-peak response and the lower the prominence of softening behavior.It can still be noted that the results reflect the very well-known phenomenon of brittleness growth under uniaxial tension with the increase in compressive strength.As expected, the similarity in the load-displacement diagrams produced with the CSMF and AMS3 corresponds to a similarity in the crack patterns generated by both strategies, as illustrated in Figure 15, for a P2-C2 sample.In this figure, the cracking process is depicted by means of crack patterns corresponding to peak and post-peak loads.The localized fracture process, which involves cracking localization and crack propagation, is clearly noticed.As can be seen, the main crack did not develop at half the height of the specimen, which is a natural consequence of the random assignment of tensile strength values.If the rupture region is not previously directed, the main crack can appear anywhere along the height of the specimen, such as what occurs in experimental tests.A natural characteristic of the explicit cracking modeling approach is that the geometrical representation of cracks depends on the degree of refinement of the mesh.Therefore, the more refined the mesh, the more cracking information can be obtained.The degree of mesh refinement adopted in the present work is in accordance with the level of cracking information required for the investigation.

Assessment of Time Savings
The efficiency of AMS3, in terms of computational execution time, was evaluated by means of discrepancy between the average simulation times related to this strategy and to the CSMF.This discrepancy is denoted by ∆ 4 .Negative values of ∆ 4 refer to time savings, and, naturally, the higher the modulus of these values, the higher the efficiency of AMS3.
A computer with the basic characteristics reported in Table 11 was used to measure the simulation times.Table 12 displays the simulated classes, together with basic information about the performed analyses: the parameters b and c, whose values, for samples related to P1 and P2, are the same as those adopted in Section 5.1; the number of Monte Carlo samples (M); and the simulation stopping criterion (SSC) employed for each class, exposed according to definitions expressed in Section 3.3.
tests.A natural characteristic of the explicit cracking modeling approach is that the geometrical representation of cracks depends on the degree of refinement of the mesh.Therefore, the more refined the mesh, the more cracking information can be obtained.The degree of mesh refinement adopted in the present work is in accordance with the level of cracking information required for the investigation.

Assessment of Time Savings
The efficiency of AMS3, in terms of computational execution time, was evaluated by means of discrepancy between the average simulation times related to this strategy and to the CSMF.This discrepancy is denoted by ∆ .Negative values of ∆ refer to time savings, and, naturally, the higher the modulus of these values, the higher the efficiency of AMS3.
A computer with the basic characteristics reported in Table 11 was used to measure the simulation times.Table 12 displays the simulated classes, together with basic information about the performed analyses: the parameters  and , whose values, for samples related to P1 and P2, are the same as those adopted in Section 5.1; the number of Monte Carlo samples (); and the simulation stopping criterion (SSC) employed for each class, exposed according to definitions expressed in Section 3.3.13 reveals the results referring to the efficiency of AMS3, which are the following: the average time required to simulate a Monte Carlo sample of the considered class, for the CSMF and AMS3; the discrepancy ∆ 4 ; and the parameter η, defined as the average percentage ratio between the number of interface elements present in the mesh at the end of the simulation with AMS3 and the maximum possible number of these elements in the analyzed case.In Table 13, the negative values of ∆ 4 indicate an efficiency of AMS3 for all the classes.The smallest moduli of ∆ 4 , 3.2% and 8.6%, are associated with very high values of η, while the largest moduli, 72.9% and 76.9%, are related to low values of η.There is coherence in this relationship between ∆ 4 and η, since high values of η denote that the amount of interface elements used in AMS3 become similar to that used in the CSMF, with AMS3 having time-consumption conditions similar to those observed in the CSMF.In this way, when the number of interface elements in AMS3 quickly approaches the maximum possible value, which happens with P1-C1 and P2-C1 (but mainly with P1-C1), this strategy shows low efficiency, tending to cause relatively small reductions in the simulation times.On the other hand, for values of η smaller than 85%, the average time savings are very significant.Thus, AMS3 proved to be noticeably advantageous, demonstrating an efficiency that naturally depends on the characteristics of the analyzed problem.
Individual results can be examined in Figure 16, which reveals the individual simulation times of 30 samples for P2-C1, P2-C2, and P2-C3, as well as the simulation times of the 20 B1-C4 samples considered in the analyses, for the CSMF and AMS3.In this figure, the difference in the simulation times for the same class of Monte Carlo samples is a natural consequence of the adopted probabilistic approach.Depending on the distribution of random tensile strength values, the crack pattern of the Monte Carlo samples can be very different, leading to distinct load-displacement diagrams (see Figure 13), similar to what happens in real experiments.Different load-displacement diagrams can be related to very different numbers of nonlinear iterations and, consequently, very different simulation times.
As can be seen in Figure 16b for samples 16 and 17, it is possible that, very infrequently, the simulation time for AMS3 is longer than that for the CSMF.The explanation for this is that, as previously justified, the stress calculation referring to the interfacial regions of solid elements does not happen in the exact same way in these strategies, so occasional differences between the mechanical behaviors produced by them can occur, mainly in the post-peak branch (see Figure 10b).Thus, the crack propagation process may not be exactly the same for both strategies, directly influencing the simulation times associated with them.For the same reason, the time reductions generated by AMS3 for some samples were considerably greater than the average reductions indicated in Table 13, as noticed, for example, in Figure 16a for sample 24 and in Figure 16b for sample 4.These observations reinforce the importance of evaluating the efficiency of AMS3 through average simulation times, favoring the statistical compensation of infrequent results.
In classes where practically the entire fracture process is concentrated in a delimited region of the mesh (P1-C3, P2-C3, and B1-C4), AMS3 saved a marked amount of time for absolutely all the analyzed samples, in accordance with what can be seen in Figure 16c for samples P2-C3 and in Figure 16d for samples B1-C4.It is important to note that the concentrated fracture process had different causes: predominantly sudden failure, for samples P2-C3, for example, as illustrated in Figure 10c; and presence of a notch, for samples B1-C4.In classes where practically the entire fracture process is concentrated in a delimited region of the mesh (P1-C3, P2-C3, and B1-C4), AMS3 saved a marked amount of time for absolutely all the analyzed samples, in accordance with what can be seen in Figure 16c for samples P2-C3 and in Figure 16d for samples B1-C4.It is important to note that the concentrated fracture process had different causes: predominantly sudden failure, for samples P2-C3, for example, as illustrated in Figure 10c; and presence of a notch, for samples B1-C4.

Conclusions
This paper communicates the development of an adaptive mesh strategy to optimize the use of interface elements in a 3D probabilistic explicit cracking model for concrete.The

Conclusions
This paper communicates the development of an adaptive mesh strategy to optimize the use of interface elements in a 3D probabilistic explicit cracking model for concrete.The investigation involved the evaluation of three adaptive mesh strategies: AMS1, AMS2, and AMS3.Based on comparisons between mechanical behaviors, only AMS3 proved to be able to replace the classical strategy of mesh formation (CSMF) generally adopted in explicit cracking models.AMS1 and AMS2 are not suitable to replace the CSMF due to divergences in the redistribution of stresses during fracturing, which lead to remarkably different mechanical responses.On the other hand, the mechanical responses produced by using AMS3 tend to present a very high degree of equivalence with those obtained through the CSMF, as there is similarity in the redistribution of stresses.Thus, an important conclusion is that different adaptive mesh strategies can lead to distinct stress redistribution processes and, consequently, to very different mechanical responses.
Possible differences between the mechanical responses related to AMS3 and the CSMF are basically concentrated in the post-peak behavior.Such differences are generally of little importance and stem from the fact that the calculation of stresses in interfacial regions of solid elements does not happen in the exact same way for the two strategies.Nevertheless, with the application of the Monte Carlo method, the divergent behaviors tend to compensate for each other, so that there is convergence on very similar average results.
The adaptive probabilistic model based on AMS3, as well as the classical alternative, is capable of predicting the scale effect at a level similar to that experimentally observed, besides covering the occurrence of different levels of softening behavior.
Regarding the efficiency of AMS3, in terms of the computational execution time, it can be concluded that this adaptive procedure is advantageous.This strategy reduced the average simulation time for all the analyzed cases.Understandably, the level of efficiency depends on the characteristics of the considered problem.The smallest average time savings, in percentage terms, were observed for cases in which the number of interface elements quickly reached values very close to the maximum possible number.In cases where this did not occur, there were significant time savings, corresponding to reductions in the average simulation time that ranged from 32% to 77%.
It is also important to highlight that, as the probabilistic nature of the cracking model used here did not affect the development of the adaptivity strategies, AMS3 can also be applied to purely deterministic cracking models that employ interface elements.

) in which 𝑎 = 6 . 5
MPa,  = 0.35, and the values of  and  are defined by

Figure 2 .
Figure 2. Important features of the probabilistic cracking model implemented by Mota et al. [11]: use of tetrahedral solid elements that remain undamaged; use of interface elements following a brittle behavior to represent cracking; local random distribution of tensile strength values; and Monte Carlo method for structural responses.

Figure 2 .
Figure 2. Important features of the probabilistic cracking model implemented by Mota et al. [11]: use of tetrahedral solid elements that remain undamaged; use of interface elements following a brittle behavior to represent cracking; local random distribution of tensile strength values; and Monte Carlo method for structural responses.

Figure 3 .
Figure 3. (a) Mesh without interface elements with a crack detected in an interfacial region of solid elements; (b) classical strategy of mesh formation for probabilistic explicit cracking models-CSMF; (c) Adaptive Mesh Strategy 1-AMS1; (d) AMS2; (e) AMS3.

Figure 3 .
Figure 3. (a) Mesh without interface elements with a crack detected in an interfacial region of solid elements; (b) classical strategy of mesh formation for probabilistic explicit cracking models-CSMF; (c) Adaptive Mesh Strategy 1-AMS1; (d) AMS2; (e) AMS3.

Materials 2024 , 31 Figure 5 .
Figure 5. Generic load-displacement diagram, highlighting information that can be used in the simulation stopping criteria:  ,  , and  .

Figure 5 .
Figure 5. Generic load-displacement diagram, highlighting information that can be used in the simulation stopping criteria: F max , F u , and δ tot .

Figure 6 .
Figure 6.Geometries and meshes for tensile test simulation.Figure 6. Geometries and meshes for tensile test simulation.

Figure 6 .
Figure 6.Geometries and meshes for tensile test simulation.Figure 6. Geometries and meshes for tensile test simulation.

Figure 6 .
Figure 6.Geometries and meshes for tensile test simulation.

Figure 8 .
Figure 8. Mesh of the considered beam (B1).Interface elements can only be used in the region highlighted in red, whose length is 100 mm.

Figure 6 .
Figure 6.Geometries and meshes for tensile test simulation.

Figure 8 .
Figure 8. Mesh of the considered beam (B1).Interface elements can only be used in the region highlighted in red, whose length is 100 mm.

Figure 8 .
Figure 8. Mesh of the considered beam (B1).Interface elements can only be used in the region highlighted in red, whose length is 100 mm.

Materials 2024 , 31 Figure 9 .
Figure 9. Final step of the inverse analysis procedure for obtaining the most suitable pair of parameters  and .

Figure 9 .
Figure 9. Final step of the inverse analysis procedure for obtaining the most suitable pair of parameters b and c.
indicates the three interfacial regions of solid elements considered in such analyses.Regions 565 and 951 are located in the central plane of the cube, and region 1318 is inclined over region 565.Regions 565 and 1318 are located on faces of the same solid element, which corresponds to the top solid element of the interface element associated with region 565.
indicates the three interfacial regions of solid elements considered in such analyses.Regions 565 and 951 are located in the central plane of the cube, and region 1318 is inclined over region 565.Regions 565 and 1318 are located on faces of the same solid element, which corresponds to the top solid element of the interface element associated with region 565.

Figure 11 .
Figure 11.Interfacial regions of solid elements of P3 (highlighted in blue) considered in the stress redistribution verification.

Figure 11 .
Figure 11.Interfacial regions of solid elements of P3 (highlighted in blue) considered in the stress redistribution verification.

7 ,
x FOR PEER REVIEW 20 of 31 set of Monte Carlo sample curves that has characteristics that makes it close to an average curve.In the mentioned figures, the resemblance between the responses produced with the CSMF and AMS3 is clear, and occasional divergences between individual behaviors are easily compensated by the application of the Monte Carlo method.It is important to observe that parameter  exerts a strong influence on the results.Low values of this parameter are associated with high variability in the obtained diagrams.Moreover, the higher the value of , the lower the level of nonlinearity in the pre-peak response and the lower the prominence of softening behavior.It can still be noted that the results reflect the very well-known phenomenon of brittleness growth under uniaxial tension with the increase in compressive strength.

Figure 15 .
Figure 15.Crack patterns of sample P2-C2 obtained by using (a) the CSMF and (b) AMS3 for (i) the peak load and (ii,iii,iv) post-peak loads.

Figure 15 .
Figure 15.Crack patterns of sample P2-C2 obtained by using (a) the CSMF and (b) AMS3 for (i) the peak load and (ii,iii,iv) post-peak loads.

Materials 2024 ,
17, x FOR PEER REVIEW 24 of 31 through average simulation times, favoring the statistical compensation of infrequent results.

Table 1 .
Details relating to geometries and meshes of the considered computational samples.

Table 2 .
Mechanical properties and heterogeneity parameters of the simulated concretes.

Table 4 .
Results of the mechanical evaluation of the adaptive mesh strategies for samples related to P1.

Table 5 .
Results of the mechanical evaluation of the adaptive mesh strategies for samples related to P2.

Table 6 .
Stress redistribution results for a crack in region 565 of P3.

Table 6 .
Stress redistribution results for a crack in region 565 of P3.

Table 7 .
Weibull distribution parameters b and c obtained via inverse analysis.

Table 9 .
Verification of the accuracy of the inverse analyses carried out with samples B1-C4.

Table 11 .
Basic characteristics of the computer used to measure simulation times.

Table 12 .
Basic information about the analyses performed to evaluate the efficiency of AMS3.

Table 11 .
Basic characteristics of the computer used to measure simulation times.

Table 12 .
Basic information about the analyses performed to evaluate the efficiency of AMS3.

Table 13 .
Results referring to the efficiency of AMS3.